Speed of synchronization in complex networks of neural oscillators 
Analytic results based on Random Matrix Theory 
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We analyze the dynamics of networks of spiking neural oscillators. First, we present an exact 
linear stability theory of the synchronous state for networks of arbitrary connectivity. For general 
neuron rise functions, stability is determined by multiple operators, for which standard analysis is not 
suitable. We describe a general non-standard solution to the multi-operator problem. Subsequently, 
we derive a class of rise functions for which all stability operators become degenerate and standard 
eigenvalue analysis becomes a suitable tool. Interestingly, this class is found to consist of networks of 
leaky integrate and fire neurons. For random networks of inhibitory integrate-and-fire neurons, we 
then develop an analytical approach, based on the theory of random matrices, to precisely determine 
the eigenvalue distribution. This yields the asymptotic relaxation time for perturbations to the 
synchronous state which provides the characteristic time scale on which neurons can coordinate 
their activity in such networks. For networks with finite in-degree, i.e. finite number of presynaptic 
inputs per neuron, we find a speed limit to coordinating spiking activity: Even with arbitrarily 
strong interaction strengths neurons cannot synchronize faster than at a certain maximal speed 
determined by the typical in-degree. 



The individual units of many physical systems, 
from the planets of our solar system to the atoms 
in a solid, typically interact continuously in time 
and without significant delay. Thus at every in- 
stant of time such a unit is influenced by the cur- 
rent state of its interaction partners. Moreover, 
particles of many-body-systems are often consid- 
ered to have very simple lattice topology (as in 
a crystal) or no prescribed topology at all (as in 
an ideal gas). Many important biological systems 
are drastically different: their units are interact- 
ing by sending and receiving pulses at discrete in- 
stances of time. Furthermore, biological systems 
often exhibit significant delays in the couplings 
and very complicated topologies of their interac- 
tion networks. Examples of such systems include 
neurons, which interact by stereotyped electrical 
pulses called action potentials or spikes; crickets, 
which chirp to communicate acoustically; pop- 
ulations of firefiies that interact by short light 
pulses. The combination of pulse-coupling, de- 
lays, and complicated network topology formally 
makes the dynamical system to be investigated 
a high-dimensional, heterogeneous nonlinear hy- 
brid system with delays. Here we present an ex- 
act analysis of aspects of the dynamics of such 
networks in the case of simple one-dimensional 
nonlinear interacting units. These systems are 
simple models for the collective dynamics of re- 
current networks of spiking neurons. After briefiy 
presenting stability results for the synchronous 
state, we show how to use the theory of random 
matrices to analytically predict the eigenvalue 
distribution of stability matrices and thus find the 
speed of synchronization in terms of dynamical 



and network parameters. We find that networks 
of neural oscillators typically exhibit speed lim- 
its and cannot synchronize faster than a certain 
bound defined by the network topology. 



I. INTRODUCTION 

Most neurons in the human central nervous sys- 
tems communicate by sending and receiving brief stereo- 
typed electrical pulses, called action potentials or spikes. 
Via chemical synaptic connections, these spikes induce 
changes in the potential across the membrane of the 
connected postsynaptic neurons Due to this mode 
of communication, these neurons interact at discrete in- 
stances in time only - and thus behave substantially dif- 
ferent from the interacting units of many physical sys- 
tems. Other important characteristics of neuronal com- 
munication are delayed interactions (due to finite prop- 
agation speed of the spikes along axons, non-zero time 
needed for chemical processes across the synapses and 
signal transmission along the dendrites) and a complex 
wiring diagram. As in the example of neurons, many 
networks of interacting units are not arranged in regular 
lattices. Instead, single units form an intricate network 
of connections that mediate the interactions. In addition, 
these connections are often directed, meaning that a con- 
nection from one unit to another does not imply a connec- 
tion in the reverse direction. From a dynamical systems 
perspective, these aspects - discrete interaction times, in- 
teraction delays, and non-symmetric, complicated wiring 
diagram - make the theoretical investigation of the exact 
spiking dynamics of large neural networks a challenging 
task. 

Previous research has occasionally explicitely consid- 
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ered interaction delays in analytical calculations; the 
complicated topology of neural networks, however, has 
received much less attention. As a consequence, if one 
wants to uncover the dynamics beyond numerical inves- 
tigations, one is often restricted to mean field theoretical 
arguments or focuses on globally connected networks or 
on networks of simple local topology 0, S S S IE S 

Here we follow the sirnple and very useful approach 
of Mirollo and Strogatz [l2| to represent the state of a 
one-dimensional (neural) oscillator not by its membrane 
potential, but by a phase that encodes the time to the 
next spike in the absence of any interactions. In the limit 
of infinitely fast processing of incoming signals (post- 
synaptic currents), the nonhnear interactions can then 
be treated analytically in an exact manner. Following 
some previous reports 0, IT^ that used the 

advantages of the Mirollo Strogatz idea we here 
present an analytical approach to exactly determine the 
asymptotic dynamics of spiking neural networks of com- 
plicated topology. We particularly focus on how, and how 
fast, neurons can synchronize their spikes, i.e. coordinate 
their activity in time in networks of random topology. 

The paper is organized as follows. In Section II we 
briefiy introduce model networks of pulse-coupled neu- 
ral oscillators and state the research question. We are 
interested in the stability of the synchronous state and 
its asymptotic synchronization properties. Section III 
gives the details of the derivation of nonlinear strobo- 
scopic maps of perturbed synchronous states in networks 
of arbitrary connectivity. We explain the emergence of 
piecewise analytic maps where the pieces are determined 
by the temporal spiking order of a particular perturba- 
tion. This results in a multiple operator nonlinear stabil- 
ity problem. In Section IV, we derive first order operators 
from the stroboscopic maps leading to a stability opera- 
tor with multiple piecewise linear parts. Since standard 
eigenvalue analysis is not appropriate for such multiple 
operator problems, we describe an alternative method to 
demonstrate stability in Section V. Section VI shows how 
degeneracy can be enforced, i.e. how all multiple Hnear 
operators can be made degenerate to one single stabil- 
ity matrix. It turns out that the oscillator rise functions 
that guarantee degeneracy are of integrate-and-fire type. 
For this stability problem, standard eigenvalue analysis is 
suitable. For two ensembles of random networks, we first 
study their eigenvalue distributions (Section VII), ana- 
lytically predict these distributions by measures derived 
from Random Matrix Theory (Section VIII) and compare 
the results between numerics and analytics (Section IX). 
In Section X, we discuss consequences of the eigenvalue 
distributions for the speed of synchronization of neural 
oscillators. We close in Section XI, where we summarize 
the results, discuss some of their consequences and give 
a brief outlook. 

This paper presents new aspects and detailed descrip- 
tions of the determination of the asymptotic synchroniza- 
tion time by Random Matrix Theory. Parts of the results 



on stability and speed limits to network synchronization 
have been reported in brief in references jl^] and [T^ . 
respectively. Details of the stability theory, in particular 
exact eigenvalue bounds and asymptotic stability in the 
multi-operator case, not discussed here, can be found in 
[T^ . For effects on parameter inhomogeneities, leading 
to close to synchronous patterns of spikes, we refer the 
reader to flil. 



II. MODEL OF NEURAL OSCILLATORS 

Consider a system of N neural oscillators that interact 
by sending and receiving pulses via directed connections. 
The sets Pre(i) of presynaptic oscillators having input 
to an oscillator i define the network connectivity. The 
number of inputs 

fc, :=|Pre(i)| (1) 

to every oscillator i, called in-degree in graph theory j23| 
is non-zero, ki > 1, and no further restriction on the 
network topology is imposed unless otherwise stated. 

The state of an individual oscillator j is represented 
by a phase-like variable G (— oo, 1] that increases uni- 
formly in time, 

d(f>j/dt = l. (2) 

Upon crossing a firing threshold, ^j(if) > 1, at time if an 
oscillator is instantaneously reset to zero, 0j {t^) — 0, and 
a pulse is sent. After a delay time r this pulse is received 
by all oscillators i connected to j (for which j G Pre(i)) 
and induces an instantaneous phase jump 

MiU + r)+) = {U{dp,{U + r) + £,,) 

Here, Eij < are the coupling strengths from j to i, which 
are taken to be purely inhibitory (e^ < if j G Pre(i), 
Eij = otherwise) and normalized, 

N 

(3) 

i=i 

throughout this paper. 

The rise function U, which mediates the interactions, is 
monotonic increasing, U' > 0, concave (down), U" < 0, 
and represents the subthreshold dynamics of individual 
oscillators. This models the dynamics of the membrane 
potential of a biological neuron that is driven by a cur- 
rent. Note that the function U need to be defined on 
the entire range of accessible phase values. In particu- 
lar, inhibitory coupHng can lead to negative phase values 

<|>^ < 0. 

Large sparsely connected networks of inhibitory neu- 
rons were known before to exhibit irregular asynchronous 
spiking states in which excitatory drive and inhibitory 
feedback balance out and fiuctuation induce spikes ja. 
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Figure 1: Irregular, balanced activity coexists with regu- 
lar, synchronous activity. This enables switching by external 
stimulus signals. Random network of connection probabil- 
ity p = 0.2 {N = 400, I = 4.0, e = 16.0, r = 0.14). Firing 
times of five oscillators are shown in a time window At — 240. 
Vertical dashed lines mark external perturbations: (i) large 
excitatory pulses lead to synchronous state, (ii) a small ran- 
dom perturbation (|A(^i| < 0.18) is restored (iii) a sufficiently 
large random perturbation (|A<j!>i| < 0.36) leads to an irregu- 
lar state. Bottom: Time evolution of the spread of the spike 
times after perturbation (ii), total length At = 0.25 each. De- 
creasing width of the distribution indicates resynchronization. 



1^, Q|. However, in a previous study we found that 
regular states, in the homogeneous case defined by ex- 
act spike synchrony, coexist with irregular states in these 
networks at the same parameters (Fig. QJ. This means 
that by external perturbations one can switch between 
regular and irregular activity. In particular, strong exci- 
tatory synchronous inputs can synchronize the network 
activity. Strong random inputs can switch the network 
back to the balanced state. If random inputs to the syn- 
chronous state are not too strong, the activity relaxes 
back to the synchronous state. Two major questions in- 
trigued us: 1) Why, given an irregular topology of the 
network, can the regular synchronous state be stable such 
that neurons resynchronize their spikes upon sufficiently 
small perturbations? 2) What is the typical time scale 
for re-synchronization, i.e. how fast can neurons coor- 
dinate their spiking activity if they are not directly in- 
terconnected but interact on large networks of complex 
topology? 

We address these questions analytically in the follow- 
ing, focusing on the speed of synchronization. All results 
are derived for the simplest of all regular states, the syn- 
chronous periodic state, in which all neural oscillators 
exhibit identical dynamics. However, a similar approach 
can be used for cluster states in which two or more groups 
of synchronized oscillators exist 0,0], as well as for any 
periodic solution because they can be tracked analyti- 
cally in the model system used. In the presence of inho- 
mogeneity, the approach needs to be modified but similar 
principles are expected to apply. 



III. NONLINEAR STROBOSCOPIC MAPS: 
EMERGENCE OF MULTIPLE OPERATORS 

The synchronous state 

Mt)^Mt) foralH, (4) 

in which all oscillators display identical phases 0o(O on 
a periodic orbit such that 4>o{t -\- T) — (poit), is one of 
the simplest states a network of neural oscillators may 
assume. The normaHzation of the coupHng strengths (pj 
ensures that it exists but does not tell whether or not it 
is stable and an attractor of the system. To uncover this, 
we perform a stability analysis of the synchronous state 
the period of which is given by 

T r + 1 - a (5) 

where 

a^U-\U{T)+e). (6) 

For inhibitory coupling (e < 0) and sufficiently small 
delay r < 1 the total input is subthreshold, U{t) + e < 1 
such that a < 1. A perturbation 

S{0)=:S = {Si,...,Sn) (7) 
to the phases is defined by 

6,= MO) -MO)- (8) 

If we assume that the perturbation is small, in the sense 
that 

max 5i — min 5i < t (9) 

i i 

it can be considered to affect the phases of the oscillators 
at some time just after all signals have been received, 
i.e. after a time t > to + t if all oscillators have fired 
at t — to- Such a perturbation will affect the time of 
the next firing events because the larger the perturbed 
phase of an oscillator is, the earlier this oscillator reaches 
threshold and sends a signal. 

To construct a stroboscopic period-T map, 6 is ordered 
according to the rank order rank(5) of the St: For each 
oscillator i we label the perturbations Sj of its presynaptic 
oscillators j G Pre(i) according to their size 

A,^i > A,^2 > . . . > A,,fc, (10) 

where ki is the number of its presynaptic oscillators 1^ . 
The index n S {!,..., h} counts the signals that ar- 
rive successively. Thus, ii jn = jn{i) G Pre(i) labels the 
presynaptic oscillator from which i receives its n'^ signal 
during the period considered, we have 

Ai,„ = (5j,^(i) . (11) 

In addition, we define 



Ai,o = <5i ■ 



(12) 
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Figure 2: Two signals arriving almost simultaneously induce different phase changes, depending on their rank order. The figure 
illustrates a simple case where Pre(i) = and Si = 0, (a)-(c) for 6ji > Sj and (d)-(f) for Sj > Sji . (a), (d) Local patch 

of the network displaying the reception times of signals that are received by oscillator i. Whereas in (a) the signal from j' 
arrives before the signal of j, the situation in (d) is reversed, (b), (e) Identical coupling strengths induce identical jumps of 
the potential U but (c),(f) the phase jumps these signals induce are different and depend on the order of the incoming signals. 
For small \5i\ <g 1, individual phase jumps are encoded by the pi,n , see 11711 . The Figure displays an example for inhibitory 
(negative, phase-retarding) coupling but the mechanism generating multiple operators does not depend on the signs of the 
coupling strengths. 



For illustration, let us consider an oscillator i that has 
exactly two presynaptic oscillators j and / such that 
Pre(i) = {j,j'} and fci = 2 (Fig. EKjd) . For certain per- 
turbations, oscillator i first receives a signal from oscil- 
lator j' and slightly later from oscillator j. This deter- 
mines the rank order, Sj' > 5j, and hence Ai^i — Sj' and 
^i,2 = Sj (Fig.Et)- Perturbations with the opposite rank 
order, Sj > Sji, lead to the opposite labeling, A.^ i — Sj 
and Ai^2 = Sj' (Fig. [21). In general, relabeling cannot be 
achieved by permuting the indices because one oscillator 
j' may receive an input connection from yet another one 
m whereas oscillator j may not receive this connection. 

We now consider a fixed arbitrary perturbation, the 
rank order of which determines the „ according to the 
inequalities ljl()|l . Using the phase shift function ft.((/), e) — 
U~^{U{(j)) + e) and denoting 

Di^n '■— Ai^„_i — Ai^n (13) 

for n G {l,...,ki} we calculate the time evolution of 
phase-perturbations Si satisfying the bound lO , starting 
near </>o(0) = t/2 without loss of generality. The strobo- 
scopic time-T map of the perturbations. Si i— > Si{T), is 
obtained from the scheme given in Table U The time to 
threshold of oscillator i, which is given in the lower left 
entry of the scheme. 



t 







i + <5« =: 1 + A,,o 




h{T + A,i,£iji) =■■ /3i,i 




fe(/3i,l + A,2, Eija) =■ Pi,2 




ft(/3i,fcj-l + Di^ki, £ijk-) ~'- Pi,ki 




reset: 1 



Table I: Time evolution of oscillator i in response to ki suc- 
cessively incoming signals from its presynaptic oscillators jn , 
n G {1, . . . ,ki\, from which i receives the n**^ signal during 
this period. The right column gives the phases (j}i{t) at times 
t given in the left column. The time evolution is shown for 
a part of one period ranging from (f>i ~ t/2 to reset, 1-^0, 
such that (j>i — in the last row. The first row gives the ini- 
tial condition (j}i{0) = t/2 + Si . The following rows describe 
the reception of the h signals during this period whereby the 
phases are mapped to /9i,n after the n"^ si gnal has been re- 
ceived. The last row describes the reset at threshold such that 
the respective time T^'^'' — t/2 — Ai_fc. -|- 1 — Pi,ki gives the 
time to threshold of oscillator i. 



is about (/'o(O) = t/2 smaller than the period T. Hence 
the period-T map of the perturbation can be expressed 
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as 



[T) = T- T^'> --= -a + A,.,, (15) 
where a is given by Eq. 0. 



IV. MULTIPLE FIRST ORDER OPERATORS 

In order to perform a local stability analysis, we con- 
sider the first order approximations of the maps de- 
rived in the previous section. Expan ding fii^ki for small 
Di^n ^ 1 one can proof by induction [T^l that 



ki 



(16) 



n=l 



where 



U'{U-\U{T)+e)) ^ > 

for n e {0, 1, . . . , ki] encodes the effect of an individual 
incoming signal of strength £ij„. The statement x = 
y means that x = y + J2i,n^i^i,n) A,n 0. 

Substituting the first order approximation Eq. l[T6|l into 
Eq. lO using Eq. leads to 



S^{T) ^ ^p,,„_i(A»,„_i - A^,„) + A,,fc, (18) 



such that after rewriting 



(19) 



ra=l 



i,n — f^i^i) for n G 



to first order in all Ai „. Since A 
{1, . . . ,ki} and A^^o = Si according to Eqs. l(lT|l and ltT2ll . 
this results in a first order map 

5{T) = AS (20) 

where the elements of the matrix A are given by 

Pi,n - Vi,n-i if j = in G Pre(i) 
P»,o if j = i (21) 

ifj ^ Pre(i) U {i}. 

As for the nonlinear stroboscopic maps Ijl5|l . because j„ 
in Eq. Il2ip identifies the v}^ pulse received during this 
period by oscillator i, the first order operator depends 
on the rank order of the perturbations, A = A(rank(J)). 
The variables pi_„ encode phase jumps evoked by all 
pulses up to the nth one received. Since the matrix ele- 
ments lf2T|l are differences of these Pi^n , matrix elements 
Aij and Aijr with j ^ j' have in general different values 
depending on the order of incoming signals. 



This multi-operator problem is induced by the struc- 
ture of the network together with the pulsed interactions. 
For networks with homogeneous, global coupling different 
matrices A can be identified by an appropriate permuta- 
tion of the oscillator indices. In general, however, this 
is impossible. Thus even for a network of given number 
of neuronal oscillators at given connection strengths and 
given delay and interaction function, the stability of the 
synchronous state is described by many different opera- 
tors that depend on the rank order of the perturbation. 



V. ALTERNATIVE METHOD 
TO DETERMINE STABILITY 

In most stability problems for periodic orbits in dy- 
namical systems theory, finding the eigenvalues of an ap- 
propriate stroboscopic map is sufficient for determining 
the stability of the orbit. Typically, one eigenvalue equals 
one and corresponds to perturbation along the periodic 
orbit trajectory such that there is no restoring force. If all 
other eigenvalues are smaller than one in absolute value, 
the periodic orbit is asymptotically stable and all suffi- 
ciently close initial states converge to it. 

On the contrary, the multi-operator property of the 
stability problem considered here implies that standard 
eigenvalue analysis fails. However, we found other meth- 
ods to determine the stability of the synchronous periodic 
state. We present the results briefiy in the following. 

To show plain (non-asymptotic) linear stability, ob- 
serve that the row-sums of the stability matrices are nor- 
malized, 



N 



(22) 



refiecting the invariance of the periodic orbit with re- 
spect to perturbations along it. Given that the coupling 
strengths are purely inhibitory, < , one can show 
that the pi_„ (Eq. IjlTIl l are positive and bounded above 
by one. 



< V^,n < 1, 
and that they increase with n, 

Pi^n — l ^ Pi^n • 



(23) 



(24) 



Hence the nonzero off-diagonal elements are positive, 
= Pi^n - Pi,n-i > such that 



> 



(25) 



for alH, j G {1, . . . , N}. Moreover the diagonal elements 



A, 



Pz,0 



U'{t) 



U'{U-^U{T)+e)) 



An 



(26) 
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are identical for all i and satisfy 



< An < 1 



(27) 



because U is monotonically increasing, C/'(0) > 0, and 
concave down, U" {(j)) < 0, for all (f>. It is important to 
note that A has the properties Eqs. lf22|l - ll27jl indepen- 
dent of the parameters, the network connectivity, and 
the specific perturbation considered. With these observa- 
tions, it is straightforward to show that the synchronous 
state is stable in the sense that small perturbations can- 
not grow: To first order, a given perturbation S = ^(0) 
satisfies 



max 



|(5(r)|| := max|(S,(T) 

N 

max V I Ay II I 

j 

max \Aij I max \ 5k 
j 

max Aij max | Sk \ 



max \Sk\ 



< 



< 



where we use the vector norm 



max 1 5i 



(28) 
(29) 
(30) 
(31) 

(32) 

(33) 
(34) 

(35) 



Thus the length of a perturbation vector does not in- 
crease during one period implying that it does not in- 
crease for an arbitrary long time. Using methods from 
graph theory, one can show that for strongly con- 
nected networks (in which every oscillator can be reached 
from every other by following a directed path on the net- 
work) the synchronous state is asymptotically stable such 
that from sufficiently close initial conditions the spiking 
activity will become exactly synchronous. The results 
on asymptotic stability use the recurrence properties of 
strongly connected networks and rely on the fact that ev- 
ery oscillator can communicate with every other at least 
indirectly. Thus the results for plain and asymptotic sta- 
bility are independent of the specific choice of parame- 
ters, eij < 0, T € (0,1), the potential function U{(f>), 
and the rank order of the perturbation. They are de- 
rived without using the eigenvalues or eigenvectors of a 
given stability matrix and solve the stability problem ex- 
actly. In summary this means that any network of the 
type described above, with normalized inhibitory cou- 
pling Eq. ^ exhibits a synchronous state that is at least 
marginally stable; it is moreover asymptotically stable if 
the network is strongly connected. 



A simple intuitive argument why networks of inhibito- 
rily coupled neural oscillators synchronize can be ob- 
tained from the response dynamics of individual units. 
Fig. El If two (or more) neurons simultaneously receive 
inhibitory input of the same size, their potential is de- 
creased by the same amount such that their potential 
difference stays unchanged. Due to the negative curva- 
ture of the rise function that mediates the negative input, 
this, however, leads to a decrease of their phase differ- 
ences, which encode the future spike times. This intu- 
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Figure 3: Intuitive synchronization mechanism: Inhibition 
synchronizes due to the concavity of U. Simultaneously re- 
ceived inhibitory input decreases phase differences between 
the receiving oscillators, {Mi^) - </'i(*'^)| > IM^^) - M't)\ ■ 



itive explanation holds for simple situations like globally 
coupled systems with homogeneous coupling strengths. 
However, the synchronization dynamics is more compli- 
cated if the inputs are not of equal size or only one input 
exists for some unit, as e.g., in a ring of neurons. 

In the case of integrate and fire rise functions U = Uip 
one can derive [;l3| stability results based on the eigen- 
system because all stability operators become degener- 
ate; see below for details of the degeneracy for networks 
of integrate and fire neurons. The Gersgorin disk the- 
orem then bounds all eigenvalues in a disk of radius 
tq = 1 — Aq centered at Aq, touching the unit circle from 
the inside at z = 1. It ensures that the eigenvalue largest 
in magnitude is Ai = 1, with corresponding eigenvector 
t>i cx (1, 1, . . . , 1)^. For strongly connected networks, the 
Perron Frobenius theorem implies that this eigenvalue is 
unique, i. e. all other eigenvalues are smaller than one 
in absolute value. Thus all perturbations that contain 
components other than Vi will decay towards a uniform 
perturbation. Such an analysis confirms that networks 
of arbitrary connectivity are at least marginally stable 
and strongly connected networks exhibit asymptotically 
stable synchronous states, as shown above by alternative 
methods. 

If networks consist of several strongly connected com- 
ponents, the analysis is much more involved and struc- 
tural identification of strongly connected components and 
the wiring among them is required. Such networks dis- 
play a novel kind of non-synchronous activity that is con- 
trolled by the coarse and fine scale structure of the net- 
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work, cf. |2j. This state seems to be universal among 
networks of coupled oscillators exhibiting a synchroniza- 
tion mechanism. 



VI. ENFORCING DEGENERACY: 
THE PHOENIX INTEGRATE- AND-FIRE 

Prom the general class of concave increasing rise func- 
tions, we now derive a subclass of rise functions in which 
all multiple operators degenerate to a single stability ma- 
trix if the coupling strength are suitably chosen. Inter- 
estingly, it turns out that the class of standard leaky 
integrate-and-fire oscillators provides potential functions 
consistent with this condition. 

A general potential function U that is monotonically 
increasing, U'[4') > 0, and concave (down), U"{(j)) < 0, 
yielded stability operators A in the first order map Ij2()|l 
that are defined by their respective matrix elements Ij21|l 
in terms of differences of the pi,„ (Eq. ITTjl that in turn 
describe the effect of the nth signal received by oscillator 
i within the period considered. Thus, the actual stability 
operator to be used for a specific perturbation depends 
on the rank order of the incoming signals given this per- 
turbation. Can the multiple linear operators be made 
degenerate? If so, the eigensystem of the resulting ma- 
trix completely described the asymptotic synchronization 
dynamics. 

Consider a network for which the coupling strengths of 
all presynaptic oscillators j £ Pre(j) are identical. 



(36) 



for each oscillator i. For such a network, two matrix ele- 
ments are interchanged at the boundary of the domains of 
definition of an individual operator. For instance, assume 
that an oscillator i has exactly two presynaptic oscillators 
j and j'. If a perturbation is changed such that Sj > Sj' 
is turned into Sj < 5ji , the operator A will change from 
A = A'^^^ to A — A^^) where the non-zero off-diagonal 
elements of row i read 



Af^ ^ Aij^ -Pifi 



A. 



(2) 



A,, 



Pi,2 - PiS 



Pd37) 

p^,m 



respectively. As above, ji labels the oscillator presynap- 
tic to i that has sent the first signal to i during the period 
considered, and j2 labels the presynaptic oscillator that 
has sent the second one such that 

ji = j and j2 = j' ^ 5j > Sj> , (39) 
ji = / and j2 ^ 6j' > Sj . (40) 

Degeneracy of these two, in general distinct, cases re- 
quires that 



A. 



(fc) 4 ^(0 



(41) 



for k,l £ {1,2} or, equivalently. 



Pi,2 -PiS =Pt,l -Ptfl ■ 



(42) 



If every oscillator i £ {1, . . . , N} in the network has ki 
presynaptic oscillators, the degeneracy condition is eas- 
ily generalized to Eq. Ij41|l with k and I running over 
all different stability matrices that occur for all possible 
differently ordered perturbations. Expressed in terms of 
the pi^n, which describe the effect of individual incoming 
pulses, we obtain 



Pi. 



' Pi,n—1 — Pi,'i 



' Pi,m — 1 



(43) 



for alH £ {1, . . . , N} and all n, m £ {1, . . . , ki}. 
If we define 



U'iU-^UiT)+e)) 



where 



E 



m—1 



ne 

ki 



(45) 



for n < ki, the requirement l(i3|l is satisfied if 

q'{x) = const (46) 

in the relevant interval x £ [e, 0]. Note that e < because 
we consider inhibitory coupHng. The first derivative of 
q{x) satisfies 



U"iU-'(U{T)+x)) ^ U"{h{x)) 
^ ^ ' U'(U-^{U{t)+x)) • U'{h{x)) 



(47) 



where h{x) = U~^{U{t) + x) is an invertible function 
of X. Together with Eq. Ij46|l this leads to a differential 
equation 



U" = cU' 



(48) 



where c £ M is a constant. The solution [/((/)) — a + be'^'^ 
with constants a, 6, c £ M together with the normalization 
C/(0) = 0, U{1) — 1, and the monotonicity and concavity 
requirements, U'{4>) > and U" {(p) < 0, yield the one- 
parameter family of solutions in integrate-and-fire form 

C/(0) = (7iF(0)-/(l-e-^^'^) (49) 

where / > 1 and Tip = ln(//(/ - 1)) > 0. This leads to 

C/i'f(</)) ^/Tipe-^^-, (50) 



f/fp^(y) = ^ln(l-f)"\ 



(51) 



and 
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such that 

UiAUrFiUicj)) + e)) = Tip (/e"*^- - e) (53) 

and 



Pi,n 



/g-rTiF _ £ 



(54) 
(55) 



Thus, by construction, if we substitute Sij^ = e/kt all 
non-zero off-diagonal elements 



1 e 



/g-TTiF - s ki 

in one row i of the stability matrix are identical, 



A - — A- 



(56) 



(57) 



for all n, m £ {1, . . . , ki}. 

One should note that, given the coupling strengths sat- 
isfy Eq. Ij36ll ■ the condition Ij46|l is sufficient but not nec- 
essary for degeneracy of all operators. At given parame- 
ters and a given network connectivity, one can construct 
potential functions that fulfill condition l(i6ll only on (lo- 
cal) average such that the requirement for identical (non- 
zero) off-diagonal matrix elements in each row Ij43|l is still 
satisfied. If we do not a priori fix the parameters and the 
network structure, however, the potential function Uif 
uniquely leads to operator degeneracy within the class of 
concave down, increasing functions. 

This degeneracy has important consequences: 
Whereas for the multi-operator problem the dynamics 
in the vicinity of the synchronous state is determined 
by an (unknown but deterministic) sequence of different 
linear operators, the dynamics in case of degeneracy 
is determined by the eigenvectors and eigenvalues of 
a single matrix A. In particular, the second largest 
eigenvalue 



Am := max{|Ai| : |Ai| < 1} 



(58) 



of this matrix A determines the asymptotic speed of con- 
vergence towards the synchronous state. 



mn + l)T)\r.Xl\SilT)\ 



(59) 



VII. LOCATION OF EIGENVALUES 
IN LARGE RANDOM NETWORKS 



How fast do random networks synchronize? The char- 
acteristic asymptotic time of synchronization, r^yn — 
— l/ln(Am), see Eq. Ij93ll below, is given in terms of the 
second largest eigenvalue Am that we determine from the 
distribution of eigenvalues in the following sections. In 
this section, we present examples for the distribution of 
eigenvalues of stability matrices describing the asymp- 
totic dynamics of large asymmetric random networks of 
integrate-and-fire oscillators in the vicinity of the syn- 
chronous state. From the details of the analysis de- 
scribed above, we know that all eigenvalues must be lo- 
cated in a Gersgorin disk K in the complex plane that 
is centered at < 1 (Eq- I27|l and has radius 1 — 
such that it contacts the unit circle at z = 1 from the 
inside. In the following, we consider neural oscillators 
that interact inhibitorily on two classes of random net- 
works (defined in subsections IVII Al and IVTl B|l . The po- 
tential functions of the oscillators are of the integrate- 
and-fire form U{<j)) = UiF{(f>) = /(I - e--^^'"), where 
Tip — ln(//(/ — 1)). The non-zero coupling strength 
are chosen according to = e/ki. We consider only 
sparsely connected networks which lead to sparse stabil- 
ity matrices where we term a matrix "sparse" if at least a 
positive fraction of its entries is zero in the limit of large 
N. 



A. Networks with constant in-degree 

The first class of networks is given by random networks 
in which all oscillators i have the same number ki = k 
of presynaptic oscillators which are independently drawn 
from the set of all other oscillators with uniform proba- 
bility. When increasing the network size N, the number 
of connections k per oscillator is kept fixed. We numeri- 
cally determined the eigenvalues of different stability ma- 
trices changing the network size N S {2^, . . . , 2^^}, the 
in-degree fee {2, . . . 2^}, and the dynamical parameters 
e, T, and / such that Aq € [0.6,0.9]. In general, we 
find that, for sufficiently large N and sufficiently large k, 
the non-trivial eigenvalues resemble a disk in the complex 
plane that is centered at about but has a radius r that 
is smaller than the upper bound given by the Gersgorin 
theorem 



for n,l ^ 1. 

Interestingly, the derivation of a condition for degener- 
acy led to the standard leaky integrate-and-fire model as 
a subclass of models that imply degeneracy for suitably 
chosen coupling strengths. Starting from this degenerate 
case of operators now enables us to develop a charac- 
terization of the synchronization dynamics in terms of 
eigenvalues of that operator. 



r<l-Ao. (60) 

Note that, due to the invariance of the periodic orbit with 
respect to globally constant phase shifts, there is always a 
trivial eigenvalue Ai = 1. As an example, the eigenvalue 
distributions in the complex plane are displayed in Fig. 31 
for specific parameters and differently sized networks. 
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Figure 4: Distribution of eigenvalues Xi in the complex plane for networks of fixed in-degree fc = 8 and different sizes (a) 
N = 32, (b) A'' = 128, (c) A'" = 512. For large networks, the non-trivial eigenvalues seem to be distributed uniformly on a disk 
in the complex plane. The arc through the trivial eigenvalue (dot at 2: = Ai = 1) is a sector of the unit circle. Parameters of 
integrate-and-fire oscillators are / = 1.1, e = —0.2, r = 0.05. 
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Figure 5: Distribution of eigenvalues Xi in the complex plane for networks of fixed connection probability p = 0.1 and different 
sizes (a) A'' = 32, (b) A'' = 128, (c) A'' = 512. For large networks, the non-trivial eigenvalues seem to be distributed uniformly 
on a disk in the complex plane, the radius of which shrinks with increasing network size. The arc through the trivial eigenvalue 
(dot at 2; = Ai = 1) is a sector of the unit circle. Parameters of integrate-and-fire oscillators are / — 1.1, e = —0.2, r = 0.05. 



B. Networks with constant connection probability VIII. PREDICTIONS FROM RANDOM 

MATRIX THEORY 



The second class of networks is given by random net- 
works for which every connection between any oscillator 
i and any other oscillator j ^ i is present with given 
probability p. When increasing N, this probability is 
kept fixed such that the number of connections per os- 
cillator is proportional to N. As for the other class of 
random networks, we find numerically that the distri- 
bution of non-trivial eigenvalues resemble disks in the 
complex plane that are smaller than the Gersgorin disk 
iffiTl but centered at about the same point ^o- We nu- 
merically determined the distribution of eigenvalues for 
N e {2^ . . . , 214}, p g [0.01, 0.2], and the dynamical pa- 
rameters £, T, and / such that ^0 S [0.6,0.9]. Figure 
displays examples of eigenvalue distributions for differ- 
ently sized networks at the otherwise identical parame- 
ters. 



The results of the previous section indicate that the 
eigenvalues of stability matrices for large asymmetric ran- 
dom networks of integrate-and-fire oscillators are located 
in disks in the complex plane if the network size N is 
sufficiently large. If this could be demonstrated indepen- 
dent of specific parameters, it would be guaranteed that 
all non-trivial eigenvalues are separated from the unit cir- 
cle. Thus the main condition required for the robustness 
of the stable synchronous state under a structural pertur- 
bation to the dynamics of the system would be satisfied. 
Moreover, the asymptotic synchronization time can be 
predicted analytically from these results. 

How can we predict the location of the eigenvalues? 
Since we are considering random networks, a natural 
starting point is the theory of random matrices. Random 
Matrix Theory has been investigated intensively since 
the early 1950s [Hi (see also [IH, Ull) and turned out 
to be a valuable tool for both qualitative and quanti- 
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tative description of spectral properties of complex sys- 
tems. For instance, it describes level correlations in nu- 
clear physics Jj^ as well as quantum mechanical aspects 
of chaos [26, 27|. In applications of Random Matrix The- 
ory to physical problems, it is generally assumed that 
the details of the physical system are less important for 
many statistical properties of interest. Often it turns out 
that important statistical properties such as the distribu- 
tion of spacings of energy levels in quantum systems are 
well described by the respective properties of random ma- 
trices that respect the same symmetries as the physical 
system. Both theoretical investigations and applications 
of Random Matrix Theory have focused on symmetric 
matrices. Asymmetric matrices are less well understood 
and found only limited applicability. In the following, we 
will evaluate the applicability of Random Matrix Theory 
for estimating distributions of eigenvalues of asymmetric 
stability matrices. 



Ensembles of symmetric and asymmetric 
random matrices 



For the case of real symmetric random TV x A^-matrices 
J — with independent, identically distributed compo- 
nents Jij — Jji , it is believed |2M l2il | that there are ex- 
actly two universality classes. Every ensemble of matrices 
within one of these universality classes exhibits the same 
distribution of eigenvalues in the limit of large matrices, 

— > oo, but the eigenvalue distributions are in general 
different for the two classes. Both universality classes 
are characterized by specific ensembles of matrices the 
elements of which are distributed according to a simple 
probability distribution. The class of sparse matrices is 
represented by the probability distribution 



for A" ^ 1. For the Gaussian symmetric ensemble, it 
is known (23, I3| that the distribution of eigenvalues 
/'Gauss ('^) i'^ the limit A^ ^ cx) is given by Wigner's semi- 
circle law 



PGauss('^) 



2^(4r2-A2)i if |A|<2r 







otherwise. 



(65) 



The ensemble of sparse matrices [23,[23,I13,E3 exhibits a 
different eigenvalue distribution Psparse('^) that depends 
on the finite number k of nonzero entries per row and 
approaches the distribution PGauss('^) the limit of large 
k such that 



lim pi 

k — >QO 



sparse 



(A) 

P Gauss 

(A). 



(66) 



It is important to note that in the limit of large A^ the 
distributions Pspai-ge ^^'^ Pcauss eigenvalues depend only 
on the one parameter r, that is derived from the variance 
of the matrix elements Ij64|l . 

For real, asymmetric matrices (independent Jij and 
Jji), there are no analytical results for the case of sparse 
matrices but only for the case of Gaussian random ma- 
trices. The Gaussian asymmetric ensemble (e.g. Eq. Ij62|l 
with independent Jij and Jji) yields the distribution 
of complex eigenvalues in a disk in the complex plane 



PGauss(A) — 




if |A| < r 
otherwise 



(67) 



where r from Eq. II64|I is the radius of the disk that is 
centered at zero. Like in the case of symmetric matrices, 
this distribution also depends only on one parameter r, 
that is derived from the variance of the matrix elements. 



k f 1 

Pspa.rse{Jij) jy*^ ( ^ 



(61) 



B. Stability matrices and the Gaussian asymmetric 
ensemble 



where k is the (finite) average number of entries in any 
row i and S{-) is the Dirac delta distribution. The class of 
Gaussian random matrices is exemplified by a Gaussian 
distribution of matrix elements 



PGa 



^(Jj^) = N^{2tts'^)~^ exp 



NJl 
2s2 



(62) 



To obtain symmetric matrices, one chooses Jij = Jji and 
Jii = for both ensembles. Thus the arithmetic mean of 
the eigenvalues is zero, 

N N 

[H:=^E^«-]^E^'^ = (63) 



and the ensemble variance of the matrix elements scale 
like 



(64) 



In the numerical studies of stability matrices for ran- 
dom networks (Sec. lVIlfl . we observed that all non-trivial 
eigenvalues of sparse stability matrices A are located on 
or near a disk in the complex plane (Figures ^ and . 
Since this is also predicted by the theory of asymmetric 
Gaussian random matrices, let us compare these predic- 
tions to numerical results. If the distribution of eigen- 
values of sparse asymmetric random matrices Pgparse fo'^ 
/c 3> 1 is approximately equal to the distribution of Gaus- 
sian asymmetric matrices, Psparse(A) ~ PGauss(A)) in anal- 
ogy to the case of symmetric matrices Ij66|l , and Random 
Matrix Theory is applicable to the stability matrices at 
all, we can obtain an analytical prediction for the radii 
of the disks of eigenvalues. 

The elements of the original stability matrix A have 
an average 



1 

TV 



(68) 
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and a second moment 



i=i 



V 



(69) 



where the off-diagonal sum is bounded above and below 
by 



AT 



<Y.Ai<{i^Aor 



in £1X7 kj 



due to the normalization Ij22|l . 

The variance cr^ — [Afj] — [Aij]^ given by 



(70) 



TV 



E 



42 



N 



1 

7V2 



is thus also bounded 



N N{N - 1) 7V2 



1 



(71) 



1 

iV2 

(72) 



because max^ ki < N — 1. The eigenvalues of the original 
matrix A have the average value 



1 ^ 1 ^ 



(73) 



To directly compare the ensemble of the stability matri- 
ces considered here to random matrices with zero average 
eigenvalue, (A^) — 0, and given variance H64|l . we trans- 
form the stability matrix A to 



A',i = A,.j - AqS, 



(74) 



for I G {1, . . . , N}. Here 6ij denotes the Kronecker delta, 
Sij = 1 if i = j and 5ij = if z 7^ j. The transformation 
to A' shifts all eigenvalues by —Aq and hence the average 
value of the eigenvalues to 



In addition 



and 



Ao Ao) 



N 



/|2 

N 



such that the variance is 



4 



1 

N 



42 
N 


2Ao 

iV2 


/|2 
^0 

iV2 


N 

E4 


(1 


TV 



A 



(75) 
(76) 

(77) 

(78) 
(79) 



The eigenvalue distribution of this ensemble of rescaled 
stability matrices A' for random networks may be com- 
pared to the Gaussian asymmetric ensemble with zero 
average eigenvalue and variance a\i . In such a compar- 
ison, the additional eigenvalue Ai = 1 of A, is neglected. 
This should not matter for large networks {N 3> 1). 

It is important to note that we compare the location 
of eigenvalues of a sparse matrix with deterministic non- 
zero entries at certain random positions with the eigen- 
value distribution of the Gaussian ensemble, which con- 
sists of fully occupied matrices with purely random en- 
tries. 

If we assume that the eigenvalue distributions for these 
two ensembles of networks with fixed in-degree and net- 
works with a fixed connection probability are similar to 
those for random matrices, we obtain a prediction 



(80) 



for the radius of the disk of eigenvalues from Eq. H64II . 
For further investigations, we consider the two exem- 
plary classes of large random networks of integrate-and- 
fire oscillators discussed in Sec. IVIII If we assume 
that the stability matrix A has exactly k non-zero off- 
diagonal elements per row and identical coupling strength 
Eij = e/k between the integrate-and-fire oscillators, the 
off-diagonal sum is exactly equal to 



E4.- 



such that the variance of A equals 

N 



Nk 



and the variance of A! is given by 



4. 



1 

772 



1 

iV 



(81) 



(82) 



(83) 



If we now take the prediction from Random Matrix The- 
ory rRMT for the radius r of the disk of eigenvalues of the 
stability matrices, we obtain 



''RMT 



N^-OA' = (1 - Ao) 



1 

TV 



(84) 



In random networks where all oscillators have exactly k 
presynaptic oscillators, the approximation for the vari- 
ance of A (and thus of A!) is exact. If the random net- 
work is constructed by choosing every connection inde- 
pendently with probability p, the variance l(82|l is only an 
approximation because we replaced by k~^ which 

gives the order of magnitude of the number of connec- 
tions as a function of N . 

Substituting ^0 = 1 — Ej j^i fo'" integrate and fire 
neurons (j56ll into the radius prediction trmt, Eq- II84|I . 
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we obtain 

which explicitely contains all parameters of the system. 

IX. NUMERICAL TESTS OF EIGENVALUE 
PREDICTIONS 

We verified this scaling law for different parameters 
Aq determined by different /, e, and r and found good 
agreement with numerically determined eigenvalue dis- 
tributions. We compared the theoretical prediction Ij84|l 
to the numerical data for both ensembles considered in 
Sec.rvTfl 

A. Networks with constant in-degree 

At a given network connectivity and given parame- 
ters, we obtained all eigenvalues of the stability matrix 
A for several network sizes N and in-degrees k. We find 
that the prediction obtained from Random Matrix The- 
ory well describes the numerically determined eigenval- 
ues. Examples of eigenvalue distributions for matrices at 
fixed k and three different N are shown in Fig. |2l 

There are several ways to numerically estimate the ra- 
dius of the disk of eigenvalues. For illustration, we use 
three different estimators here. The real part estimator 

''Re := 7: ( maxRe(Ai) - minRe(Ai) ) (86) 

estimates the radius from the maximum spread of eigen- 
values parallel to the real axis. Typically, r^e should give 
an estimate that is too low compared to the radius ob- 
tained from the eigenvalues of an ensemble of matrices 
because it measures the maximal spread in one direction 
only. This is circumvented by the radial estimator 

r^ad := max I A, - {Ao - (1 - Ao)N~^)\ (87) 
1^1 

that finds the maximum distance of any non-trivial eigen- 
values from the average of the non-trivial eigenvalues, 
(A,)^_^i = Ao - (1- Aa)N-'^ + 0{N-^). This estima- 
tor should yield an approximation that may be too large 
compared to the respective ensemble average. The aver- 
age estimator 

S 1 ^ 

r.. := - {Ao - {1 - Ao)N-^)\ (88) 

i=2 

estimates the radius r of a circle from the average dis- 
tance (d) of eigenvalues from its center, because 

{d) = / r'^p{r')drd(p = -r (89) 

Jo Jo 3 



if we assume a uniform p{r') = l/(7rr^) for r' < r and 
p(r') = otherwise llfiTji . This estimate has the advan- 
tage, that it contains information from all eigenvalues in 
contradistinction to the two other estimators. Its dis- 
advantage is that one has to assume a priori a uniform 
distribution of non-trivial eigenvalues. As displayed in 
Fig. all three estimators converge towards the radius 
predicted by the random matrix model for large N and 
given in-degree k. Varying the in-degree k at fixed TV also 
yields excellent agreement between the numerical data 
and the theoretical predictions for sufficiently large N 
and k. An example is displayed in Fig. |H1 

For both, networks of fixed k and networks of fixed 
p, there are deviations for small and even for intermedi- 
ate N, because the prediction trmt was obtained from 
Random Matrix Theory that is exact only in the limit 
N ^ (X, and the finite-size scaling of trmt was assumed 
to resemble the scaling of the variance of finite matrices. 
Furthermore, as discussed above, the numerical estima- 
tors of the radius rely on assumptions that are fulfilled 
only approximately. For sufficiently large networks, how- 
ever, the theoretical prediction agrees well with the nu- 
merical data. 

Thus there is a gap of size 

g=l-Ao-rao (90) 

between the non-trivial eigenvalues for large networks 
and the unit circle, where 

Too lim TRMT = (1 - Ao)k-^/^. (91) 

This indicates that the stability of the synchronous state 
in the model system considered is robust, i.e., sufficiently 
small perturbations to the systems dynamics will not al- 
ter the stability results. 

Nevertheless, there is an important restriction to these 
results. Given a fixed in-degree k, the limit N co 
is not described by the theory derived in the previous 
section, because the structure of the network considered 
and thus the structure of the stability matrices is only 
well defined if the network is not connected in the sense 
that every oscillator has at least one presynaptic oscilla- 
tor. However, the probability that at least one oscillator 
is disconnected from the remaining network approaches 
one with increasing network size. Thus eigenvalue pre- 
dictions of stability matrices of networks with fixed in- 
degree k are only reasonably described for network sizes 
that are large, ^ 1, but not in the limit TV ^ 00. 

B. Networks with constant connection probability 

If we assume that every connection is present with a 
constant probability p, the network will be connected 
with probability one in the limit N ^ 00 because the 
number of presynaptic oscillators ki follows a binomial 
distribution with average pN and standard deviation 
(p(l — p)NY/'^. In this limit, the radius of the disk of 
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Figure 6: Distribution of eigenvalues in the complex plane for networks with fixed in-degree k = 8 for different network sizes 
(a) A'' = 128, (b) A'' = 512, (c) A'^ — 4096. The disks are centered at Aq and have radius trmt, the prediction obtained 
from Random Matrix Theory. The arc through the trivial eigenvalue z = Ai = 1 is a sector of the unit circle. Parameters of 
integrate-and-fire oscillators are I = 1.1, e = —0.2, r = 0.05. 
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Figure 7: Scaling of the radius r of the disk of non-trivial 
eigenvalues with the network size A'^ at fixed in-degree fc = 32 
(1 = 1.1, e = —0.2, T = 0.05). Main panel displays the 
radius r as a function of network size A'. Symbols display 
frad (*), fav (x) and rRe (O)- Inset displays roo - r as a 
function of A^ on a doubly logarithmic scale, where Voo — 
(1 — Ao)k~^^^ . Dots display numerical data of rav. In the 
main panel and the inset, lines are the theoretical prediction 
rRMT = (l-Ao)(l/fc-l/Af)i/2. 



eigenvalues decreases with increasing network size N, see 
Fig. El 

In order to verify the scahng behavior of the radius of 
the eigenvalue disk for large stability matrices A, we nu- 
merically determined the eigenvalues Am = max{|Ai| : 
I Ail < 1} , see Eq. lf58jl . that are second largest in abso- 
lute value. For sufHciently large N, the theoretical pre- 
diction Am ~ Aq + trmt agrees well with the numerical 
data (Fig. The radius approaches zero for large net- 
works such that the eigenvalue second largest in absolute 
value converges towards the center Aq of the disk. In 
conclusion, for large networks, all non-trivial eigenvalues 
are located near Aq and are thus bounded away from the 
unit circle. This implies that the speed of synchroniza- 
tion that is determined by Am increases with increasing 



Figure 8: Scaling of the radius r of the disk of non-trivial 
eigenvalues with the in-degree k for random networks of A^ = 
1024 oscillators (7 = 1.1, e = -0.2, t = 0.05). Main panel dis- 
plays the radius r as a function of in-degree k. Inset displays 
the same data on a doubly logarithmic scale. Symbols display 
numerical results, using the average estimator rav, lines are 
the theoretical prediction trmt = (1 — Ao)(l/k — 1/A')^''^. 

network size. Moreover, the condition necessary for ro- 
bustness against structural perturbations of the systems 
dynamics is satisfied. 



X. SYNCHRONIZATION SPEED 
AND SPEED LIMIT 

The existence of bounds on the radius of the eigenvalue 
distribution has severe consequences for the synchroniza- 
tion speed of networks of neural oscillators. Whereas 
the largest (trivial) eigenvalue Ai = 1 corresponds to 
the invariant nature of the synchronized periodic orbit, 
the second largest eigenvalue Am (Eq. EHl determines 
the asymptotic speed of synchronization starting from 
sufficiently close- by initial conditions. Because the dy- 
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Figure 9: Distribution of eigenvalues in the complex plane for networks with fixed connection probability p — 0.1 for different 
network sizes (a) A'' = 128, (b) A'^ = 512, (c) A*' = 4096. The disks are centered at Ao and have radius trmt, the prediction 
obtained from Random Matrix Theory. Note that the disk of non-trivial eigenvalues shrinks towards the point Aq in the limit 
N ^ oo. The arc through the trivial eigenvalue 2: = Ai = 1 is a sector of the unit circle. Parameters of integrate-and-fire 
oscillators are / = 1.1, e = —0.2, r — 0.05. 
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Figure 10: Maximum non-trivial eigenvalue and the radius 
of the eigenvalue distribution for random networks (same pa- 
rameters as in Fig.O. Main panel displays the maximal non- 
trivial eigenvalue Am ~ ^o + r as a function of network size A^. 
The maximal non-trivial eigenvalue converges to Aq ~ 0.83 as 
N ^ 00. Inset displays the radius r of the disk of eigenvalues 
as a function of A''. Dots display numerical results based on 
= ''rad (Eq. (|H3l), lines are the theoretical predictions for 
both, the radius r and the maximal non-trivial eigenvalue Am. 



namics can locally be approximated by a linear map, 
the synchronization of spike times is an exponential. 
Thus, denoting S'{t) :— 6{t) — lims-+oo ^{s), the distance 
A(n) := maxi \dl{nT)\/ max.; |(5^(0)| from the invariant 
state behaves as 



A(n) ~ exp(— n/rg. 



(92) 



defining a synchronization time r^yn in units of the collec- 
tive period T. The speed of synchronization t^"^ strongly 
depends on the parameters. For instance, as might be ex- 
pected, synchronization is faster for stronger coupling. 

Given the results from Random Matrix Theory derived 
above, we can deduce an expression for the synchroniza- 



'syn 



-l/ln(A„) 
-l/ln(^o +?'rmt) 



(93) 



from the prediction of the second largest eigenvalue ~ 
^0 + ''RMT- In general, upon increasing the coupHng 
strength e, the center Aq of the disks of eigenvalues is 
moved towards zero, as can be seen from the defining 
equation Ij26|l . This means that, as expected, stronger in- 
teraction strengths lead to faster synchronization of the 
spiking activity if all other dynamical and network pa- 
rameters are kept fixed. There is, however, a speed limit 
to synchronization if the in-degree of the network is fi- 
nite. The speed limit plays a noticeable role if the typical 
in-degree k is significantly smaller than the number TV of 
oscillators in the network. 

For networks with constant in-degree (Fig. Illll , the 
radius of the eigenvalue disk converges to a positive con- 
stant with increasing network size N. This means that 
the second largest eigenvalue does not converge to zero 
as the coupling increases arbitrarily strong. Thus the 
asymptotic synchronization time ll93)l is bounded below 

by 



syn 



2 

Ink 



iVln(fc) 



(94) 



for large TV and thus limited by the network connectivity 
(cf. the asymptotes in Fig. [TT)l . 

For networks with fixed connection probability p, the 
radius of the eigenvalue disk does converge to zero with 
increasing network size TV. However, for any finite net- 
work with finite number of connections per oscillator, it 
has a positive radius and again leads to a speed limit, 
see Fig. Note that in the example displayed the typ- 
ical number of connections per oscillator is as large as 
k « pN « 409 but the speed Hmit is still prevalent. 

Can we intuitively understand the speed limit that is 
enforced by the topology of the network, parameterized 
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Figure 11: Asymptotic synchronization time in random net- 
works, (a) Network with fixed in-degree h = k = 32 
{N = 1024, I = 1.1, T = 0.05, e^j = e/k for j G Post(i)) . 
The inset shows the distance A of a perturbation 5 from the 
synchronous state versus the number of periods n {s = —0.4). 
Its slope yields the synchronization time Tsyn shown in the 
main panel as a function of coupling strength |e|. Simulation 

data (O): theoretical prediction ( ) derived in this paper, 

its infinite coupling strength asymptote ( ). (b) Network 

of A'' = 2048 neural oscillators and connection probability 
p = 0.2 . Other parameters and inset as in (a). Note that 
although the typical in-degree is changed drastically from (a) 
to (b), the synchronization speed limit is hardly affected. 



by its typical in-degree? Consider a large number of neu- 
ral oscillators connected via a network of complicated 
topology. If from the fully synchronous state (Fig. \T% ) 
only one oscillator is perturbed away (Fig. \V2b ) this 
constitutes a simple example of resynchronization. One 
might imagine that all the other oscillators are pulling 
the phase of the perturbed one back to their common 
phase fFig. 112b). This would explain why, with increas- 
ing coupling strengths, synchronization would be faster 
- the stronger the local puUing force, the faster the lo- 
cal resynchronization. If that was the only mechanism 
involved, the network could be resynchronized arbitrar- 
ily fast using sufficiently large coupling strengths. The 
actual mechanism, however, is non-local. Because in the 
linearized dynamics each neural oscillator performs lo- 
cal averaging, see Eqs. (I21ll - lj27|l . of their own phase and 



Figure 12: Schematic illustrating the mechanism of resyn- 
chronization in a network of pulse-coupled neural oscillators. 
A collection of oscillators (connections not shown) at specific 
phases illustrated as "time on the clocks", (a) Unperturbed, 
fully synchronous state, (b) One oscillator perturbed (out of 
phase), (c) Purely local restoring of phases might seem to 
be the natural way for resynchronization but it is not possi- 
ble because local averaging of phases implies spreading of the 
perturbation such that finally (d) all oscillators of the network 
agree on a common phase that does not equal their original 
one. 



those phases of its presynaptic oscillators, the common 
phase of the resynchronized state will be globally agreed 
on fFig.ll2H). i.e. determined by the phases of all oscilla- 
tors in the network. Neural oscillators can only interact 
with their neighbors, and, due to their pulsed interac- 
tions, only at discrete times once a period. For inhibitory 
interactions this means that the time between communi- 
cation events is bounded below by Atinteract > inde- 
pendent of the delay time r. At long times, the averaging 
has to be performed all over the network, thus restricting 
the speed of synchronization. 



XI. CONCLUSION 

We have investigated the dynamics of synchronization 
in networks of neural oscillators with complex connec- 
tion topology. We first described the stability analysis 
for the general case and found that the arising nonlinear 
and first order mappings have multiple state dependent 
parts. As an important consequence standard eigenvalue 
analysis of the first order system is not suitable. Us- 
ing alternative methods, we demonstrated that the sim- 
plest periodic state, the synchronous state in which all 
neurons fire periodically at identical times, is stable for 
inhibitory coupHng, independent of the specific network 
topology. Second, to study the speed of synchronization, 
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we derived a subclass of models for which all parts of the 
first order stability operators become degenerate. This 
class in general requires rise functions of integrate-and- 
fire type. Subsequently, we used Random Matrix Theory 
to analytically predict the speed of synchronization via 
the eigenvalue distributions depending on dynamical and 
network parameters. Numerical estimates are in excellent 
agreement with our theoretical predictions. 

Although the theory used is based on Gaussian 
(i.e. fully occupied) matrices in the limit TV ^ oo, our 
results also hold for sparse random networks with mod- 
erately large finite N. Moreover, it is known that the 
eigenvalue distribution of the sparse symmetric random 
matrix ensemble converges towards the eigenvalue distri- 
bution of the Gaussian symmetric ensemble in the limit 
/c — > oo. It is not clear whether a similar relation holds 
for Gaussian and sparse asymmetric ensembles as well. 
In fact it is an open question why the the Gaussian en- 
semble actually describes the synchronization of sparse 
random networks even for small k « 10^ rather than 
only for fc — > cx). 

Our results also indicate that stable synchrony is com- 
mon to a class of neural oscillators and not restricted to 
the specific model considered here. Moreover, given the 
expression for the speed of synchronization, we discov- 
ered a speed limit to synchronization on networks that is 
controlled by the typical in-degree of each oscillator, i.e. 
the number of other oscillators it receives input from. 
The dependence of the speed limit on the in-degree is 
logarithmic such that even for large in-degree the speed 
limit is significant. 

The application of Random Matrix Theory in the 
present study suggests that it might well be possible to 
analytically predict dynamical properties of other sys- 
tems from their structure, using an ansatz comparable 
to ours. Examples of the application of Random Ma- 
trix Theory in ecology are provided in references 



They were restricted to the dynamics near fixed points. 
Due to the idealization in the model class considered here, 
it was possible to analytically predict dynamical aspects 
near invariant (periodic) solutions that are not simple 
fixed points using Random Matrix Theory. 

Some straightforward generalizations of possible appli- 
cation include less simple periodic states like cluster peri- 
odic orbits [4, 10, 17] or periodic patterns of spikes which 
occur in the presence of heterogeneity |^. 
More interesting, and certainly more involved possibili- 
ties for Random Matrix Theory applications may arise 
if the dynamics becomes unstable. Of particular interest 
for theoretical neuroscience may be saddle periodic orbits 
which imply a high degree of flexibility when switching 
between states [lO, 'sg^'i^, IH El Ea]. Starting from 
the class of systems considered in the current paper, the 
next step into this direction would be to consider orbits 
that arise in networks where inhibitory and excitatory 
recurrent interactions coexist [5, 6, 7]. 

Our approach is not restricted to the well known Er- 
dos Renyi random graphs considered here. If other net- 
work topologies have to be considered, we expect that 
under some additional assumptions, just the associated 
random matrix ensemble could be used to describe the 
linearized dynamics of such systems. For future investi- 
gations of synchronization properties of networks, scale 
free and small world | l44l45i | topologies constitute promis- 
ing candidates because these networks might be analyti- 
cally tractable but nevertheless appear to reflect impor- 
tant aspects of real world networks. 

We thank Michael Denker, Tsampikos Kottos, Peter 
Miiller, Haim Sompolinsky, Martin Weigt and Annette 
Zippelius for useful discussions and comments on this 
work. Supported in part by the Federal Ministry of Ed- 
ucation and Research (BMBF), Germany, under grant 
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Neuroscience BCCN Gottingen, Project A2). 



[1] F. Rieke, D. Warland and W. Bialek, Spikes: Exploring 
the Neural Code (MIT Press, Cambridge, Massachusetts, 
1999). 

[2] P. Bressloff, S. Coombes and B. de Souza, Phys. Rev. 

Lett. 79, 2791 (1997). 
[3] C. van Vreeswijk, L. F. Abbott and G. B. Ermentrout, 

J. Comp. Neurosci. 1, 313 (1994). 
[4] U. Ernst, K. Pawelzik and T. Geisel, Phys. Rev. Lett. 

74, 1570 (1995). 
[5] C. van Vreeswijk and H. Sompolinsky, Science 274, 1724 

(1996). 

[6] C. van Vreeswijk and H. Sompolinsky, Neural Comput. 

10, 1321 (1998). 
[7] N. Brunei and V. Hakim, Neural Comput. 11, 1621 

(1999). 

[8] C. van Vreeswijk, Phys. Rev. Lett. 84, 5110 (2000). 
[9] D. Hansel and G. Mato, Phys. Rev. Lett. 86, 4175 (2001). 
[10] M. Timme, F. Wolf and T. Geisel, Phys. Rev. Lett. 89, 
154105 (2002). 



[11] N. B. A. Roxin and D. Hansel, Phys. Rev. Lett. 94, 
238103 (2005). 

[12] R. E. MiroUo and S. H. Strogatz, SIAM J. Appl. Math. 

50, 1645 (1990). 
[13] U. Ernst, K. Pawelzik and T. Geisel, Phys. Rev. E 57, 

2150 (1998). 

[14] W. Senn and R. Urbanczik, SIAM J. Appl. Math. 61, 
1143 (2001). 

[15] M. Timme, F. Wolf and T. Geisel, Phys. Rev. Lett. 89, 
258701 (2002). 

[16] M. Timme, F. Wolf and T. Geisel, Phys. Rev. Lett. 93, 
074101 (2004). 

[17] M. Timme, F. Wolf and T. Geisel, Chaos 13, 377 (2003). 

[18] M. Timme and F. Wolf, in preparation (2005). 

[19] M. Denker, M. Timme, M. Diesmann, F. Wolf and T. 

Geisel, Phys. Rev. Lett. 92, (2004). 
[20] G. Chartrand and L. Lesniak, Graphs and Digraphs, 

3. edition (Chapman and Hall, Boca Raton, Florida, 

1996). 



17 



[21] M. Timme, submitted (2005). 

[22] E. P. Wigner, Proc. Camb. Phil. Soc. 47, 790 (1951). 
[23] Statistical Theory of Spectra: Fluctuations, ed. C. E. 

Porter (Academic Press, New York, 1965). 
[24] M. L. Mehta, Random Matrices (Academic Press, New 

York, 1991). 

[25] T. Guhr, A. Miiller-Groeling and H. A. Weidenmiiller, 

Phys. Rep. 4-6, 189 (1998). 
[26] O. Bohigas, M. J. Giannoni and C. Schmit, Phys. Rev. 

Lett. 52, 1 (1984). 
[27] F. Haake, Quantum Signatures of Chaos (Springer, New 

York, 2001). 

[28] A. D. Mirlin and Y. V. Fyodorov, J. Phys. A 24, 2273 
(1991). 

[29] Y. V. Fydorov and A. D. Mirlin, Phys. Rev. Lett. 67, 
2049 (1991). 

[30] A. J. Bray and G. J. Rodgers, Phys. Rev. B 38, 11461 

(1988). 

[31] G. J. Rodgers and A. J. Bray, Phys. Rev. B 37, 3557 
(1988). 

[32] V. L. Girko, Theory Probab. Appl. 29, 694 (1985). 

[33] H. Sommers, A. Crisanti, H. Sompolinsky and Y. Stein, 



Phys. Rev. Lett. 60, 1895 (1988). 
[34] R. M. May, Nature 261, 459 (1976). 
[35] V. K. .Jirsa and M. Ding, Phys. Rev. Lett. 93, (2004). 
[36] D. Z. Jin, Phys. Rev. Lett. 89, 208102 (2002). 
[37] C. Borgers and N. Kopell, Neural Comput. 15, 509 

(2003). 

[38] R.-M. Memmesheimer and M. Timme, in preparation 
(2005). 

[39] D. Hansel, G. Mato and C. Meunier, Phys. Rev. E 48, 
3470 (1993). 

[40] M. Rabinovich, A. Volkovskii, P. Lecanda, R. Huerta, 
H. D. 1. Abarbanel and G. Laurent, Phys. Rev. Lett. 87, 

068102 (2001). 

[41] A. Zumdieck, M. Timme, T. Geisel and F. Wolf, Phys. 

Rev. Lett. 93, (2004). 
[42] P. Ashwin and M. Timme, Nonlinearity 18, 2035 (2005). 
[43] P. Ashwin and M. Timme, Nature 436, 36 (2005). 
[44] R. Monasson, Eur. Phys. J. B 12, 555 (1999). 
[45] I. J. Farkas, I. Dertoyi, A.-L. Barabftsi and T. Vicsek, 

Phys. Rev. E 64, 026704 (2001). 



